SS4A Data Analysis
library(tidyverse)
library(tidycensus)
library(sf)
library(gt)
library(mapview)
# Make sure to read the GEOID as a char!
# Run `get_towns_in_region.R` to get this table.
# We need this for a lot of our work, so we'll load it up early!
towns_in_MAPC <- read_csv(
"data/processed/towns_in_MAPC_BRMPO.csv",
col_types = list(col_character()))Step 1: Fatality Count
We downloaded data from the Fatal
Analysis Reporting System (FARS) FTP site for the years 2016 through
2020. We then filtered a specific table within those files
(accident.csv) for those crashes that occurred in
Massachusetts. We combined those files together for the entire state. We
start from that simplified file.
I manually checked that Wenham, Hamilton, Norfolk, Rockport, and Essex do indeed have 0 fatal crashes during this time period. They do not show up in the accident files. Several towns were miscoded (5 = Acton, 170 = Cambridge, 125 = Braintree, 1065 = Reading, and one unclassified town was Malden). )
We also have to clean up our datasets so the town names are the same. They didn’t supply town ids we can join to a census dataset. We need to join on cityname (or in this case, the altered city_join field).
acc_ma_16_20 <- read_csv("./data/processed/accidents_ma_2016_2020.csv", show_col_types = FALSE) |>
select(STATE, COUNTY, COUNTYNAME, CITY, CITYNAME, FATALS, fars_year) |>
mutate(city_join = case_when(
CITYNAME == "BOXBORO" ~ "BOXBOROUGH",
CITYNAME == "FOXBORO" ~ "FOXBOROUGH",
CITYNAME == "MANCHESTER" ~ "MANCHESTER-BY-THE-SEA",
CITYNAME == "BROOKLINE VILLAGE" ~ "BROOKLINE",
CITYNAME == "CHARLTON DEPOT" ~ "CHARLTON",
CITYNAME == "EAST BILLERICA" ~ "BILLERICA",
TRUE ~ CITYNAME
))
# Zero Fatality Munis: Wenham, Hamilton, Norfolk, Rockport, EssexWe can group the data together by municipality and sum the data to find out the 5-year statewide totals. We’ll then need to pare it back to the MPO region.
acc_bymuni <- acc_ma_16_20 |>
group_by(STATE, COUNTYNAME, city_join) |>
summarize(total_fatals_16_20 = sum(FATALS))
acc_MAPC <- acc_bymuni |>
full_join(towns_in_MAPC |>
mutate(NAME_short = toupper(NAME_short)),
by = c("city_join" = "NAME_short")) |>
filter(!is.na(GEOID)) |>
mutate(total_fatals_16_20 = replace_na(total_fatals_16_20, 0))
acc_MAPC_summ <- acc_MAPC |>
group_by(RPA, MPO) |>
summarize(total_fatals_16_20 = sum(total_fatals_16_20)) |>
mutate(annual_fatal_rate = total_fatals_16_20 / 5)
acc_MAPC_summ |>
ungroup() |>
janitor::adorn_totals() |>
gt() |>
tab_header(title = "Fatal Crashes in MAPC Communities",
subtitle = "as Reported in the 2016-2020 FARS") |>
fmt_number(columns = total_fatals_16_20, decimals = 0) |>
fmt_number(columns = annual_fatal_rate, decimals = 1) |>
cols_label(total_fatals_16_20 = "Total Fatalities (5 Years)") |>
cols_label(annual_fatal_rate = "Average Annual Fatalities") |>
tab_source_note(source_note = "Note: RPA = Regional Planning Agency. MPO = Metropolitan Planning Organization. MAPC = Metropolitan Area Plannning Council. OCPC = Old Colony Planning Council. FARS = Fatality Analysis Reporting System") |>
tab_source_note(source_note = "Source: 2016-2020 FARS Accident Tables (FATALS), Municipalities by RPA and MPO by CTPS")| Fatal Crashes in MAPC Communities | |||
|---|---|---|---|
| as Reported in the 2016-2020 FARS | |||
| RPA | MPO | Total Fatalities (5 Years) | Average Annual Fatalities |
| MAPC | Boston | 542 | 108.4 |
| MAPC/OCPC | Old Colony | 15 | 3.0 |
| Total | - | 557 | 111.4 |
| Note: RPA = Regional Planning Agency. MPO = Metropolitan Planning Organization. MAPC = Metropolitan Area Plannning Council. OCPC = Old Colony Planning Council. FARS = Fatality Analysis Reporting System | |||
| Source: 2016-2020 FARS Accident Tables (FATALS), Municipalities by RPA and MPO by CTPS | |||
Information from MassDOT’s
Impact Portal suggest that 534 fatalities occurred in the MAPC
region: 
Step 2: 5-Year Average Fatality Rate
Download 2020 US Census Data
We will use the census API to download population data from the 2020 US Census.
vars <- tidycensus::load_variables(year = 2020, dataset = "pl")
ma_towns <- get_decennial(
geography = "county subdivision",
variables = "P1_001N",
year = 2020,
state = "MA",
geometry = TRUE,
cb = FALSE
)
MAPC_towns <- ma_towns |> left_join(towns_in_MAPC, by = "GEOID") |>
filter(str_detect(RPA, "MAPC"))Let’s make a quick map to check if that looks right.
MAPC_towns |> mapview(zcol = "RPA")Summarize 2020 Census Data for the Region
pop_MAPC_summ <- MAPC_towns |>
st_drop_geometry() |>
group_by(RPA, MPO) |>
summarize(total_pop = sum(value)) |>
ungroup()
pop_MAPC_summ |>
janitor::adorn_totals() |>
gt() |>
tab_header(title = "Population of MAPC and the Boston MPO",
subtitle = "Based on the 2020 U.S. Census") |>
fmt_number(columns = total_pop, decimals = 0) |>
cols_label(total_pop = "Total Population") |>
tab_source_note(source_note = "Note: RPA = Regional Planning Agency. MPO = Metropolitan Planning Organization. MAPC = Metropolitan Area Plannning Council. OCPC = Old Colony Planning Council") |>
tab_source_note(source_note = "Source: 2020 U.S. Census: P.L. 94-171, Table: P1 (Total Population, P1_001N)")| Population of MAPC and the Boston MPO | ||
|---|---|---|
| Based on the 2020 U.S. Census | ||
| RPA | MPO | Total Population |
| MAPC | Boston | 3,357,194 |
| MAPC/OCPC | Old Colony | 78,565 |
| Total | - | 3,435,759 |
| Note: RPA = Regional Planning Agency. MPO = Metropolitan Planning Organization. MAPC = Metropolitan Area Plannning Council. OCPC = Old Colony Planning Council | ||
| Source: 2020 U.S. Census: P.L. 94-171, Table: P1 (Total Population, P1_001N) | ||
Combine Fatalities and Population
Now we can merge the fatalities and the population to estimate the fatality rate.
fatal_rate <- left_join(acc_MAPC_summ,
pop_MAPC_summ,
by = c("RPA", "MPO")) |>
janitor::adorn_totals() |>
mutate(
fatal_per_100k = annual_fatal_rate/total_pop * 100000)
fatal_rate |> gt() |>
tab_header(title = "Fatality Rate in MAPC Communities",
subtitle = "Based on the 2020 U.S. Census and 2016-2020 FARS data") |>
fmt_number(columns = total_pop, decimals = 0) |>
fmt_number(columns = total_fatals_16_20, decimals = 0) |>
fmt_number(columns = annual_fatal_rate, decimals = 1) |>
fmt_number(columns = fatal_per_100k, decimals = 2) |>
cols_label(total_fatals_16_20 = "Total Fatalities (5 Years)") |>
cols_label(annual_fatal_rate = "Average Annual Fatalities") |>
cols_label(total_pop = "Total Population") |>
cols_label(fatal_per_100k = "5-Year Average Fatality Rate (per 100,000 persons)") |>
tab_source_note(source_note = "Note: RPA = Regional Planning Agency. MPO = Metropolitan Planning Organization. MAPC = Metropolitan Area Plannning Council. OCPC = Old Colony Planning Council. FARS = Fatality Analysis Reporting System") |>
tab_source_note(source_note = "Population source: 2020 U.S. Census: P.L. 94-171, Table: P1 (Total Population, P1_001N)") |>
tab_source_note(source_note = "Fatality source: 2016-2020 FARS Accident Tables (FATALS), Municipalities by RPA and MPO by CTPS")| Fatality Rate in MAPC Communities | |||||
|---|---|---|---|---|---|
| Based on the 2020 U.S. Census and 2016-2020 FARS data | |||||
| RPA | MPO | Total Fatalities (5 Years) | Average Annual Fatalities | Total Population | 5-Year Average Fatality Rate (per 100,000 persons) |
| MAPC | Boston | 542 | 108.4 | 3,357,194 | 3.23 |
| MAPC/OCPC | Old Colony | 15 | 3.0 | 78,565 | 3.82 |
| Total | - | 557 | 111.4 | 3,435,759 | 3.24 |
| Note: RPA = Regional Planning Agency. MPO = Metropolitan Planning Organization. MAPC = Metropolitan Area Plannning Council. OCPC = Old Colony Planning Council. FARS = Fatality Analysis Reporting System | |||||
| Population source: 2020 U.S. Census: P.L. 94-171, Table: P1 (Total Population, P1_001N) | |||||
| Fatality source: 2016-2020 FARS Accident Tables (FATALS), Municipalities by RPA and MPO by CTPS | |||||
Step 3: Underserved Communities
Process Underserved Community Data
We can go to the Areas of Persistent Poverty Project (APP) and Historically Disadvantaged Community (HDC) Status Tool referenced in the NOFO to obtain underserved census tracts.
The data needs to be cleaned before it’s joined to population data. We start by using the 2010 Census geography as that matches the format of the downloaded data. We’ll identify the places that are DACs using the 2010 geography, then we’ll handle bringing that to 2020 when we have it a little more manageable. Census tracts split and merge over time as the population changes. The Census Bureau likes to keep tracts around 4,000 people.
underserved <- read_csv("./data/base/RAISE_Persistent_Poverty.csv") |>
janitor::clean_names()
underserved_filter <- underserved |>
filter(a_state == "Massachusetts") |>
mutate(underserved = e_app_county_meets_definition == "Yes" |
f_app_census_tract_meets_definition == "Yes" |
g_hdc_census_tract_meets_definition == "Yes")
# The detailed census tract data comes with similar information attached.
# We can use that information to build the join. First we need to download the full data.
# It appears that the definition file used by USDOT for the DACs was based on the 2010 Census
# geography.
ma_tracts_10 <- get_decennial(
geography = "tract",
variables = "P001001",
year = 2010,
state = "MA",
geometry = TRUE,
cb = FALSE,
keep_geo_vars = TRUE)
ma_tracts_10_mod <- ma_tracts_10 |>
separate(NAME, into = c("tract_name", "county_name", "state_name"), sep = ",") |>
mutate(tract_name = trimws(tract_name),
county_name = trimws(county_name)) |>
left_join(underserved_filter,
by = c("tract_name" = "c_census_tract_name",
"county_name" = "b_county")
)
ma_tracts_10_mod <- ma_tracts_10_mod |>
filter(county_name %in% c("Plymouth County", "Norfolk County", "Worcester County",
"Middlesex County", "Essex County", "Suffolk County")) |>
filter(underserved == TRUE)Find the Underserved Tracts in MAPC/Boston MPO
We clip the dataset to contain just those that are within the
region’s town boundaries. I was having trouble getting the intersect
feature to work as is–I suspect it has something to do with the borders
of tracts being counted as an intersection. st_covered_by
seems to work though.
underserved_inMAPC <- st_join(ma_tracts_10_mod,
MAPC_towns |> group_by(RPA, MPO) |> summarize(),
join = st_covered_by) |>
filter(str_detect(RPA, "MAPC"))
underserved_inMAPC |> select(GEOID10, e_app_county_meets_definition,
f_app_census_tract_meets_definition,
g_hdc_census_tract_meets_definition,
RPA, MPO) |> st_write("./data/processed/underserved_inMAPC_2010CensusTracts.shp",
delete_dsn = TRUE)## Warning in abbreviate_shapefile_names(obj): Field names abbreviated for ESRI
## Shapefile driver## Warning in CPL_write_ogr(obj, dsn, layer, driver,
## as.character(dataset_options), : GDAL Error 1: ./data/processed/
## underserved_inMAPC_2010CensusTracts.shp does not appear to be a file or
## directory.## Deleting source `./data/processed/underserved_inMAPC_2010CensusTracts.shp' failed
## Writing layer `underserved_inMAPC_2010CensusTracts' to data source
## `./data/processed/underserved_inMAPC_2010CensusTracts.shp' using driver `ESRI Shapefile'
## Writing 119 features with 6 fields and geometry type Multi Polygon.underserved_inMAPC |> select(GEOID10, e_app_county_meets_definition,
f_app_census_tract_meets_definition,
g_hdc_census_tract_meets_definition,
RPA, MPO) |> select(GEOID10, RPA, MPO) |> st_write("./data/processed/underserved_inMAPC_2010CensusTracts.geojson",
delete_dsn = TRUE)## Deleting source `./data/processed/underserved_inMAPC_2010CensusTracts.geojson' failed
## Writing layer `underserved_inMAPC_2010CensusTracts' to data source
## `./data/processed/underserved_inMAPC_2010CensusTracts.geojson' using driver `GeoJSON'
## Writing 119 features with 3 fields and geometry type Multi Polygon.mapview(ma_tracts_10_mod, col.regions = "black") +
mapview(underserved_inMAPC, col.regions = "red") +
mapview(MAPC_towns, zcol = "MPO")underserved_inMAPC |>
st_drop_geometry() |>
group_by(RPA, MPO) |>
summarize(DAC_total = sum(value)) |>
janitor::adorn_totals() |>
gt() |>
fmt_number(columns = DAC_total, decimals = 0) |>
cols_label(DAC_total = "DAC Population") |>
tab_header(title = "Population Living in Underserved Census Tracts in Massachusetts (2010 Census)",
subtitle = "Based on the USDOT's definition of disadvantaged communities ('DACs') from the Justice40 Initiative") |>
tab_source_note(source_note = "Note: RPA = Regional Planning Agency. MPO = Metropolitan Planning Organization. MAPC = Metropolitan Area Plannning Council. OCPC = Old Colony Planning Council. DAC = Disadvantaged Community") |>
tab_source_note(source_note = "Classification source: Areas of Persistent Poverty Project (APP) and Historically Disadvantaged Community (HDC) Status Tool: https://datahub.transportation.gov/stories/s/tsyd-k6ij") |>
tab_source_note(source_note = "Population source: 2010 US Census: SF1, P001001 (Total Population)")| Population Living in Underserved Census Tracts in Massachusetts (2010 Census) | ||
|---|---|---|
| Based on the USDOT's definition of disadvantaged communities ('DACs') from the Justice40 Initiative | ||
| RPA | MPO | DAC Population |
| MAPC | Boston | 475,609 |
| MAPC/OCPC | Old Colony | 1,753 |
| Total | - | 477,362 |
| Note: RPA = Regional Planning Agency. MPO = Metropolitan Planning Organization. MAPC = Metropolitan Area Plannning Council. OCPC = Old Colony Planning Council. DAC = Disadvantaged Community | ||
| Classification source: Areas of Persistent Poverty Project (APP) and Historically Disadvantaged Community (HDC) Status Tool: https://datahub.transportation.gov/stories/s/tsyd-k6ij | ||
| Population source: 2010 US Census: SF1, P001001 (Total Population) | ||
That’s useful to know, but it’s still the wrong year. We need to go figure out how to deal with the changes between 2010 and 2020.
ma_tracts_20 <- get_decennial(
geography = "tract",
variables = "P1_001N",
year = 2020,
state = "MA",
geometry = TRUE,
cb = FALSE,
keep_geo_vars = TRUE)
mapview(ma_tracts_20, col.regions = "grey50") + mapview(underserved_inMAPC, col.regions = "red")Hmm. So our options are to do some sort of allocation based on area
or perhaps population. We could use
tidycensus::interpolate_pw() might have what we need. Or we
can basically say, we’ll take the tracts that essentially match up with
what they were trying to do with the initial dataset. We can say any
tract that’s MOSTLY covered by one of the tracts will be
included. Then we can take a peek and see if it looks good.
# We'll start by selecting tracts that touch our tracts. It'll make things quicker.
ma_tracts_20_nearby <- st_intersection(ma_tracts_20, underserved_inMAPC) |>
st_set_geometry(NULL) |>
select(GEOID) |>
distinct() |>
pull() ## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries# No idea why tidyverse isnt working here. But we want the area of each tract.
# We'll use this to figure out which 2020 tracts to pull in.
ma_tracts_20_near <- ma_tracts_20 |>
filter(GEOID %in% ma_tracts_20_nearby)
# Find the area of the full 2020 Census tracts. This will be our total area.
ma_tracts_20_near$area_tract <- st_area(ma_tracts_20_near)
mapview(ma_tracts_20_near) + mapview(underserved_inMAPC, col.regions = "red")# Intersect the underserved areas from 2010. We're going to use this intersection
# to figure out the proportion of each new tract that's close to the underserved areas.
ma_tracts_20_filter <- st_intersection(ma_tracts_20_near,
underserved_inMAPC |> select(GEOID10 = GEOID)) |>
st_make_valid() # This seems important with some odd geometry that results from that.## Warning: attribute variables are assumed to be spatially constant throughout all
## geometriesma_tracts_20_filter$area_int <- st_area(ma_tracts_20_filter)
# Let's select a cutoff and see how we do. We call it filter2 because it's a pain
# if we messed up and have to go do that expensive intersection a second time.
# I'd rather use a little extra memory than wait around.
ma_tracts_20_filter2 <- ma_tracts_20_filter |>
group_by(GEOID) |>
mutate(pct_area = units::drop_units(area_int / area_tract)) |>
ungroup() |>
st_make_valid() |>
filter(pct_area > 0.2) |>
select(GEOID) |>
st_drop_geometry() |>
distinct() |> # No need for duplicates here!
pull()Ok! Now all that’s left to do is see how that affects our 2020 Census selection. That looks pretty good! We’re not including a sliver of Waltham, but we are including a little extra in the other direction. The left maps show the 2020 Census tracts and the right map show how that compares to the 2010 geography. This seems to keep most of the idea behind the DACs.
m1 <- mapview(ma_tracts_20, col.regions = "grey80", lwd = 1.5) +
mapview(ma_tracts_20 |> filter(GEOID %in% ma_tracts_20_filter2)) +
mapview(underserved_inMAPC, col.regions = "red")
m2 <- mapview(ma_tracts_10, col.regions = "grey80", lwd = 1.5) +
mapview(ma_tracts_20 |> filter(GEOID %in% ma_tracts_20_filter2)) +
mapview(underserved_inMAPC, col.regions = "red")
leafsync::sync(m1,m2)